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1 . Summary 


This report summarizes the results of a three year investigation 
into the structural dynamics of a cantilever turbomachine blade 
mounted on a spinning and processing rotor. 

A cumulated list of publications issued during the couTse of the 
research appear as references 1,2,3 and 5 of this final technical 
report. In addition to these references which are conference 
proceedings and journal publications, reference 17 is an interim 
report covering the first eighteen months of research under the 
Grant. The present final technical report incorporates that 
interim report, with corrections, and augments it with the results 
of the final one and one half year of effort. 

Both stability and forced vibration are considered with 
a blade model that increases in complexity (and verisindlitude) 
from a spring-restrained point mass, to a uniform cantilever, 
to a twisted uniform cantilever, to a tapered twisted cantilever 
of arbitrary cross-section. In every instance the formulation 
is from first principles using a finite element based on beam 
theory. Both ramp-type and periodic-type precessional angular 
displacements are considered. In concluding, forced vibrating 
and flutter are studied using the final and most sophisticated 
structural model. 

The analysis of stability is presented in some detail 
and a number of numerical examples are worked out. One example is 
given approximating a shroud-type restraint at the 3/4 span point. 

One other set of calculations demonstrate the role of structural 
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damping in the phenomenon. The practical occurence of this type 
of instability is discussed. 

The forced vibration problem is treated, with forcing 
present at one and two times the rotational frequency with 
amplitude dependent upon the precessional rate. Finally there 
are presented some considerations on the effect of subsonic 
aerodynamic damping on the dynamics of twisted cantilever blades. 

The conclusion is that when flutter occurs it is of the coalescent 
type with the two lowest modes coupling to account for the vanishing 
of the aerodynamic damping. 



2. Introduction to Precession-Induced Instability of Rotor Blades 

An operating turbine-type aeroengine may be considered a 
gyroscope with the turbomachine rotor blades representing a large 
number of radially disposed beams. These blades are subject to a 
wide variety of vibratory forcing mechanisms as well as flutter, 
or self-excitation. Hence rotor blade vibrations, as well as those 
of stator vanes, have received a great deal of analytical and 
experimental attention, forming a very large sub-discipline in the 
field of aeroelastic and mechanical vibrations. The excitation and 
self-excitation studied in this present research program relate to 
the fact that the engine spin axis may be forcibly processed with 
angular rates that are developed by the entire vehicle in which the 
aeroengine is mounted. 

The potential sources of precession of the rotor shaft axis 
are several in number. In highly maneuverable aircraft the pilot 
may intentionally institute rapid pull-up, nose-over or yawing 
motions with angular rates approaching, or even exceeding one radian 
per second. With a maneuver lasting only a few seconds the rotor 
may turn through hundreds of revolutions providing ample time for 
self-excitation to occur or for forced vibrations to build up. 

In addition to this ramp-type precession, the rotor may be 
subjected to an harmonic precession of the engine axis due to flying 
through a turbulent atmosphere. A wing-mounted engine would be 
subjected to the same harmonic angular rate as the wing chord in a 
flutter situation, assuming the flutter mode had an appreciable 
wing torsion component. In the latter case, however, the fluttering 
system would consist of the wing/pylon/engine and the inertial and 
aeroelastic characteristics of the engine and its nacelle would have 
to be included in the flutter stability determination. 

Another potential source of harmonic precession is the operation 
of turbine engined aircraft on rough runways, or the operation of 
land vehicles in rough terrain. In these cases, as well as in the 
turbulent atmosphere, the precession might be expected to have a 
certain statistical distribution centered about the natural frequency 
of the entire complex vibrating structure. 
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other vehicular applications of gas turbines and the seismic 
motion of stationary gas turbine mounts are additional situations in 
which gyroscopic influences on blade stability may be important. 

For the ultimate treatment of these gyroscopic phenomena in 
turboblading the aerodynamic influences should be included in the 
anaylsis. However, initial y it is convenient to ignore the aero- 
dynamic forces and deal with only elasticity and inertia. If the 
assumed operation is far from a flutter condition this assumption 
is tantamount to ignoring the aerodynamic damping provided by the 
working t^uid (air, or combustion products in a turbine). 

The order of investigation is therefore as follows: 

i) Develop a finite element for a tapered, twisted beam 

ii) Qualify the element by static loadings and rotating natural 
frequency determination 

iii) Conduct stability analyses and forced vibration analyses 

under both precession histones, steady (or ramp variation) 
and harmonic 

iv) Include the aerodynamic forces in the forced vibration 

analysis, and, if possible, in the self-excitation analysis 

v) Include the shaft restraint in the analysis (i.e., allow 

for a ’’coning" of ttie shaft as additional degrees of freedom) 

In the present report items i) to iii) are included, representing 
the results of the first 18 months of work under the present contract. 
(These results have been reported in 3 published papers 1, 2, 3). 
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3. Blade Idealized as a Cantilever Beam Under Rotor Spin and 
Space-Fixed Precession 

3. 1 Introduction 

This section deals with the study of vibratory behavior of turbo- 
machine blades under combined rotor spin and angular precession. 

The turbine blade is idealized as a cantilever beaun. Pretwist 
is not considered. The effects of blade taper, varied:>le precession 
rate, damping and Coriolis forces are studied. 

The purpose of this study is to gain insight into the com- 
plicated problem of blade vibrations under combined rotor spin and 
precession by first studying a simplified model. The results of 
this study can also be used as references for the numerical solution 
obtained from the more complete model employed in Section 4. 

The equations of motion are derived from Lagrange's equations. 

The kinetic and potential energy expressions required are developed 
in the sequel. 

3.2 Kinetic and Potential Energies 

It is assumed that the precessional velocity vector 
is space fixed and the origin of the rotor fixed xyz coordinate 
system (see Fig. 1) is a fixed inertial point. 

The angular velocity of the rotor is 


a) 


CO ( -+ _n. ( CoS 0 j -I. 


6k) 


( 1 ) 
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where and 6J is the magnitude of the rotor spin 

velocity. The position vector of a displaced point on the beam can 
be given as 

( 2 ) 

where Tx , , S? are the centroidal displacement components 

of an arbitrary beam cross section. For the purposes of this section 
Si can be considered to include only the fore-shortening effect 
(£) and is given by 

R 

Taking the time derivative of equation (2) the velocity 
of a point on the beam is 


-V 

\J = 


Sx ‘ + Sy J + S»)c ^ 


(4) 


Using equations (1) emd (2) the velocity can be written more 
explicitly 

r ^ "*■ ^x) SL sVo ® ^ j 

4 [ ^ - (x b']K 


( 5 ) 
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The kinetic energy is obtained from 

T » ^ y ^ - ^ p dv («> 

V 

We neglect all the terms which include x and y in the above 
integral because these terms are due to shear deformations and/or 
due to rotary inertia which are usually neglected in beam type 
analysis. Equation (6) can be sin^lified 

T * i ^ V • "V f 

It 

where A(z) is the variable cross sectional area. In practice 
blades usually have a very nearby linear thickness variation 


' A f (8) 

where A is the cross sectional area at the root of the blade ^ 

and t is the ratio of tip thickness to that of the base, 
r 

In order to reduce the beam problem with many degrees of 
freedom to a single degree of freedom, only bending vibration 
perpendicular to the major axis of cross section is considered. 

This leads to 




( 9 ) 



X 


Fig. 1 Space-fi: 
System 
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where ^ ie the setting angle measured from the axis of spin 
to the major axis of the local cross-section, (see Fig* 2} and 
is the displacement in the direction of the minor axis of cross 
section. Furthermore the fore- shortening effect can be written as: 



( 10 ) 


Hence the kinetic energy integral (7) can be completely expressed 
in terms of the displacement, ^ (z,t) in the flexible direction 
of the blade. For completeness we rerain the terms that contain 
cross sectional coordinates ^ and . » 

7.; . t *£• (g) * 


2 . . 




R ^ 
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( continued) 
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Neglecting shear deformations the potential energy is given as: 

u »^5 

With an known displacement function % this integral can be 
evaluated. 

3.3 Uniform Blades 
3,3.1 Equation of Motion 

To study the behavior of uniform cantilever b,i,ade8 for which 
t^ > 1, we may assume a reasonable distribution of displacements 
along the blade as follows 


- Coi 


2rJ_t-jR) 

SL 


) 


% (2, t) - a ^i) (I 


(13) 


ORIGIfJAL PiiO'i 
OF POOR QUALi iY 

u(t) being the tip displacement of the blade. This assumed form is 
inserted into the kinetic energy T and potential energy U expressions 
and the Rayleigh-Ritz procedure is applied to both integrals. 

From Lagrange's equation 



where L = T - U, the equation of the beam tip motion is obtained 


— ^ + [qo e cos ^ eo5 z © ] oT 

do 

~ s.'n e ^ 2 .U (^ ) ] 


ole 


= ~ F fci/» ^ Sind 


( 15 ) 


where the quantities ^ *2 '/ R > , ^ 

all dependent on blade parameters as follows: 


G ^ - cut 


G - < < I 


(dimensionless time) 
(typically about 0.0001) 


11,^2 

rr> L oO 


CoS p 1.^41*^ 


f.t03 — 
L 


Q, ^ -0.349 


0»7 Cos ^ 


0. %0 9 


£ 

L 


( 16 ) 
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^ Si^ ^ P 

-*■ o.S-coi’"P ■* o.«o<i-S 
_ _ I . (3<»4 Sl'/> P /L 

= Otot/L* 

F = - Losp (O.yiai + <?.«oi -2 ) 


The blade parameters are setting angle (stagger) ^ , blade 
length L, rotor radius R, mass of blade m, section modulus El. ^ 
is the irtor spin velocity whereas is the magnitude of the space 

fixed processional velocity. 

The parameters given in equation (16) are for uniform blades. 

If the blade is tapered these parameters have different expressions 
and they reflect taper effects. However, the equation of motion 
(15) retains the same form and the characteristics of the solutions 
presented in this section remain general. 

It should be pointed out here that xyz coordinates used in 

equation (15) are rotating with the rotor angular velocity tu . 

Therefore, the tip deflection u is located in a moving coordinate 

^ 1 . 

system. This in turn produces the Coriolis term ( 6 d (a ) 

which is not explicitly dependent on velocity. 

Comparing the parameter q^ to the other parameters in the 
differential equation (15) , q^ can be neglected for large blade length 
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L. Since gyroscopic inst 2 ibility is in^ortant for relatively long 
turbine or fan blades this nonlinear term will be omitted from 
further consideration in this section. For the time being it is 
worthwhile to mention that this term is due to higher order effects 
of the foreshortened axial displacement. 

The two terms on the right hand side of equation (15) represent 
forcing functions. They are generated by the centrifugal ccceleration 
associated with the rotor rotation. The first term is a time function 
with twice the frequency of the rotor spin i.e. the force becomes 
zero four times in every rotation. This force vector is in the y 
direction of the moving coordinate system. The second term is generated 
by a force vector that is in the x direction of the moving coordinate 
system. Both inhomogeneous terms on the right heind side of equation 
(15) represent components of these force vectors that cause blade 
bending in the u direction. 

3.3.2 Stability of the Linear Equation of Motion 

Equation of motion (15) is nonlinear even though the higher 
order effects of foreshortened axial displacements are ignored by 
dropping terms that are dependent on q^. In this section we shall 
omit the other nonlinear term fc ^ sin0 . Since the stability 

of the motion is determined by the homogeneous part of the equation 
of motion the linearized equation (15) is reduced to 

it — 

■ - ->r [ Oo •+ ^ OS6? -+ ^ Co3 "3 M ~ ® (17) 

The stability characteristics of this equation was discussed 


r 
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in reference (S ) . The equation is identical to the Mathieu equation 
if the terms are ignored. The stability of Mathieu 's equation 

is studied by means of the so-called Strutt diagram, (6) . The terms 
that are of order Oi^*) generate additional unstable regions and 
augment the standard stability regions in the Strutt diagram. The 
derivation of these stability curves is based on Floquet theory. 

The period of the time-dependent coefficient in equation (17) 
is 2 7r . The perturbation method of strained parameters (_7^) can be 
used to determine the periodic solutions of u. This requires the 
expansion of the solution u(0) in terms of 

U - Uo -*■ (18) 


Also it is necessary to expand the dimensionless squared 
frequency a^ in terms of . 


Qo ' bo -t et, 


(19) 


Substituting equation (18) and (19) into (17) a set of differential 
equations are obtained 


da' 


b 


O Uo 


o 


bo u , - ~ ( b>, <7 . c^s B) Uo 

de' 

~ ( b, cetB) 

etc. 


( 20 ) 
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The differential equations (20) are solved successively and 
secular terms are suppressed by imposing conditions on the expansion 
paraaneters etc. This is based on Floquet theory for the 

transition behavior between stable and unstable solutions. The 
solution of the first differential equation (20) is a harmonic function 

Oo ~ Ao Co5 ^ ^ 0 (21) 


with constants and B^. Substituting this expression into the 
second differential equation (20) 







( 22 ) 


and solving this differential equation yields 






2 


2vTb: 


_ r co5 ( > ~ ) e 

^ i 2.vfWo - 1 

^0^1 r s<o ( I — ^ Q 

2 - L "2.\J"ibD I 


CoS ( I ^ ^ 

zO^o ' 

Sib ( > ~*~ 1^0 ) ^ 
Z'Tbo ' 


] 

1 


(23) 


The first two terms on the right hand side of (23) are secular 
terms. Additional secular terms arise when Ttc - . If this 

analysis is carried out further by inserting the above solution (23) 
into the third differential equation (20) it is determined that 
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besides Jw* “ ~ ^ also gives rise to more secular 

terms. Analysis in the higher powers of ^ generates even more secular 
terms but for most practical purposes 0(€^) is accurate enough. 

Starting with - ^^7. equation (22) becomes 


c(a. 


Tjtl, 

= (-b, - 

Ao CoS a. ^ 

9 1 Ao 

39 


SL 

CoS ^ 

1 -r- 
3 ^ 


(-b, -*■ ^)Bo 


(24) 


Eliminating the secular terms in (24) requires that either 



(25) 


or 



The two values of in equations (25) and (26), when inserted 
into equation (19), correspond to two different stability transition 
curves in the plane. These curves separate stable 

from unstable regions. 

Using these results, the second order approximation of the 
transition curves can be obtained by eliminating secular terms in the 
solution of the third differential equation in (20) . This requires 
that 

( ba Q, Ao ~ O 

(27) 

and 

( ba o, <lt^) 6t> - O 
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Equations (27) and (28) indicate that 
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b 



a 


( 29 ) 


is the sufficient condition for the elimination of secular terms 
in U2* 

Equations (25), (26), and (29) are inserted into (19) to 
obtain the equations for the two transition curves that emanate from 
a^ = 1/4, e ' c 




~ ^ 






e C Qi 



(30) 


The above process is repeated at the other singular point: 

a =1 and 6^-0 . The transition curves are described by 
o 


do = I 


+ 6 


^ r lii - 

L IX 



(31) 


and 


qo I 



(32) 


Emanating from each singular point on the a^ axis of (a^,€-) 
plane there are two transition curves. The loci of transition values 
separate the (a^, ^ ) plane into regions of stability and regions of 
instability. Additional unstable regions can be derived for higher 
values of a^. However, these regions of instability are narrow 
(order of € ) and do not have much practical significance. 
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In the stable regions equation of motion (15) can be approximated 
by the following differential equation. 

— ^ *+ ZT 5* F s/rt ^ F ^ 0 

Jl9^ <33) 

This approximation is not valid in the unstable regions of the 
(a^, ) plane. 

3.3.3 Numerical Example 

The equations (30) , (31) and (32) for stability transition 
curves derived by the perturbation method contain blade parameters 
^ , R/L etc. In this section using these equations for some 

typical blade we may construct a stability chart. The stability 
chart of Fig. 12 is given for a particular set of blade parameters, 

R/L = 1.384, L = 10 in., p = -45° which yields the following para- 
meters a^ = 1, q^^ = 1, and q 2 = 2. 

For every pair of curves that emanates from the discrete 
a^ values, the region between the transition curves correspond to 
(a^, ^ ) values which give unbounded solutions of the equation of 

motion. Therefore these are regions of instability. It should be 
noted that the unstable region at 3^ = 1 is considerably narrower 
than the region at a^ = 1/4. 

We now chec)c the validity of the stability transition curves 
that are given in Fig. 12 by an independent numerical procedure. 

According to reference (£) the solutions of the homogeneous differential 
equation (17) are stable if 

I V/. 4- v/iVa-h) I < a 
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SPEED RATIO SL /u) 


Numerical results 
o Stable 
• Unstable 



Fig. 3 stability Chart for the Linear Equation of Motion 




where ( d ) and V 2 ( ^ ) are numerical solutions of (17) with 

the following initial conditions 


V, (o)-j , 

\JiCo)-^o V/a7o)-l 

j 


( 35 ) 


The predictor-corrector numerical scheme of Adams -Bash forth 
was used to obtain v^( 6 ) and V 2 ( d ) for given (a^, 6 ) values 
when a^ ® 1, = 1, and q 2 = 2, The stable and unstable points 

are mar)ced on the stability chart of Fig. 3. Excellent agreement 
is observed with the results of perturbation methods. 

3.3.4 Effect of the Nonlinear Coriolis Term 


— 

If the nonlinear Coriolis term, S'ia6 U- 

is retained in the equation of motion (15) a general stability 
analysis becomes impossible. However certain predictions can be 
made about systems which have the form 


(36) 


where 6 is small and ) is periodic with 

a period 2.7T in the product variable . In these systems 

"resonance" occurs only if 



(37) 


where p and q are small and mutually prime integers (9) . If {SZ 
is irrational and cannot be obtained by the division -of two integers 
the solutions to (36) are not periodic but usually bounded. Therefore 
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it suffices to make the stability analysis of (15) for only those 
values which can be represented as a division of two mutually 
prime integers 



( 38 ) 


Neglecting terms that depend on g^ , the equation of motion 
(15) can be written in the form of equation (36). The coefficient 
£ in this equation is the ratio of the precessional angular 
velocity SL to the spin angular velocity co . In general this 
ratio is very small compared to a^. Therefore it can be expected 
that the solution of (15) is nearly harmonic. An approximate lin- 
earized differential equation can be obtained by assuming that the 
solution consists of two parts: 


U = Ul 


( 39 ' 


The first part Uj^ satisfies the equation: 

which has the solution: 

Ul = Ai.co5>r^e® 


(40) 


(41) 


Substituting equations (39) and (40) into equation (15) and linearizing 

2 

by neglecting the ( €u) term, results in 
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-► €-c^ic^sd e'<2» ^5 a«3 

- 6 2<^a 1 M 

a ^FSio5e - <3€ 

4 SiA 0 L(l^ 

' ( € q , 9i coad -4. ^9^ Co< Qd)C4,. 

( 42 ) 


The stability of the above differential equation can be 
studied by considering the homogeneous part only. 

r Oo 4 ^ Fq, Cos^ €\z cos 2 ^ 

dfl T t 

- 6 Sto6> U * O 


( 43 ) 


Under the condition of equation (38) the period of the coefficient 
of u in the above differential equation is 2p'h' . Furthermore 

the stability of this differential equation can be studied by means 
of Floquet theory. 

Here we apply the perturbation method of strained parameters, 
again . 


M “ Uo ■*“ 6 U . -*• (44) 

Qo* bo (45) 

Substituting the above equations into (43) a set of differential 
equations are obtained 
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d^u 

de' 



*o 


j » 

— bcU» “ “(bi-*- CoS ® ■* '^^1 Sfn6 <^i^) i^e 

da* 



b. til 




“ (b, cos© - 5i;> CI^ ) (>(, 


( 46 ) 

etc. 

The differential equations (46) are solved successively 
and secular terms are suppressed by imposing conditions on the 
expansion parameters b 2 # etc. The solution of the first 
differential equation in (46) is an harmonic function 

Uo- Ac C 0 5 >Tbo ^ 'Ibc ^ ( 47 ) 

with constants and B^. Substituting this equation into the 
second differential equation in (46) 

cl'u 

^ - (b,-»-^.Co3 5 ~ -2 Si'riO Of u) Oe (46) 

da' ^ 

and solving this differential equation it is determined that secular 
terms arise when «1, = 1/2, and « 1/3. 


- 22 - 



As before, the differential equation (48) has to be solved for 
these specific values of b^ and the secular terms that arise must 
be eliminated by imposing conditions on b^. 

Starting with b^ « 1/9# equation (48) is 


U, - *■ b, Ac -I - b, 6o 


4 


3 


Ao Al 


"4 3 6«.Al 




Ao 64. 

C.S 1 



.cs-f 

1 

~ 3 

6« < 

^''>1 + 

» e e 

(49) 


Eliminating the secular terms in (49) requires that 


“b, Ao i 6 l<^3Ao ^ ^ O 

“b. 13 . ^AL^jiAo “O 


Equations (50) can be written in matrix form 



GHiCi ■ 1" ■- • ■ 
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To obtain nontrivial solutions the determinant of the coefficient 
matrix must be zero. This yields the following equation 


a 




(52) 


The two values of b^^ in equation (52) , when inserted into 
equation (45 ) , correspond to two different stability transition 
curves that emanate from a^ = 1/9 , € = 0 , 

For each b^ that is obtained from equation (52) the corresponding 
relation between and is determined from equation (51) . 

Be (53) 

When this analysis is extended to the second power of € we face 
an inconsistency that presents us from eliminating the secular terms. 
For each transition curve that corresponds to a particular b^ value 
we obtain two conditions on b 2 instead of one. 



( 54 ) 
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ta 




6 . 


- OrC/ 



( 55 ) 


Obviously conditions (54) and (55) cannot be satisfied simul- 

and B_ . Examining the two expressions for 
b^ in these equations, the difference is in the last term of the 
equations. If these terms are neglected the perturbation solution 
(44) will not be exactly periodic due to U 2 » The prerequisites for 
applying Floquet theory to the transition curve will be violated. 

The method of strained parameters fails to determine the stability 
transition curves to the second and also to higher powers of €. . 

Therefore the solution so obtained is not uniformly valid as “t -•■oo . 

Other perturbation methods have been tried without success in resolving 
this inconsistency. This remains a mathematical problem that may 
be tackled rigorously at a future time. 

Fortunately, in most practical cases the ratio, € , of 

precessional speed to rotor spin is small. The numerical error re- 
sulting from omitting these inconsistent terms is negligible. In 
order to verify these approximated stability curves the independent 
numerical method of previous section is applied to the linear equation 
(43) directly. This is presented in the next section. 

The approximate equations (54) and (55) for b 2 and equation 
(52) for bj^ are inserted into (45) to obtain the equations for the 
two transition curves that emanate from a^ = 1/9, € =0. 


taneously for arbitrary 
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(56, 57) 


The above process is repeated at other singular points: 
and a^ = 1/4, 1 etc. The approximate transition curves at a 


6=0 
= 1/4 are 



and at q = 1 are 

Qe=| -+6-^ [iBu- n/Al^ ] 



(58) 


(59) 


(60, 61) 


It should be mentioned that these are not the only stability 
curves. Additional unstable regions can be derived for higher values 
of a^. However, these additional regions are exceedingly narrow and there- 
fore they are less important. 
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All of these transition curves contain and which are 
the initial condition constants of the basic solution u , The 

Xi 

dependence of stability on the initial conditions is typical 
for nonlinear systems. But these transition curves do not depend 
on the constants and B^ in equation (47) of the perturbed solution. 
This indicates that the stability curves are valid for any arbitrary 
initially perturbed motion. 

3.3.5 Numerical Example 

The equations for stability transition curves (56) , (57) , (58) , 
(59), (60), and (61) contain the initial condition constants 
A^ and Bj^ as well as the blade parameters p , R/L, L etc. For 
each of these parameters one can construct a family of stability 
curves while the other parameters are held constant. We choose to 
examine the effect of initial conditions on the stability curves 
for a given set of blade parameters. 

Stability charts of Figs. 4 , 5 , 6 are constructed for the 

same blade data that were used in section 4.3.3 which yield a^ = 1 , 

^1 ~ ^2 ~ *^3 “ 0.1547. 

For this purpose two dimensionless initial value parameters, 
and y^-3. are defined as follows: 


is the initial displacement, u_ (0) and B_ is proportional 

1j Li 

to the initial velocity. 

Fig. 4 , Fig. 5 , and Fig. 6 correspond to different initial 
displacement values, >7, =0, =0.05, and =0.1 respectively. 

In each figure four pairs of transition curves are drawn for four 
different values of : 0, 0.05, 0.1, and 0.5. Note that =0.5 

corresponds to a significant initial speed. For example, if the spin 
velocity is 4000 rpm the initial speed for the transition curves 
emanating from a^ = 1/4 is 87 ft/sec. 

In the unstable regions, between a pair of curves emanating 
from discrete a^ values, the perturbed motion u grows and consequently, 
according to equation (39), the blade tip motion u is unstable. 

In Fig. 4 the transition curves for rj, - O are of 
particular interest. These are stability transition curves when the 
nonlinear Coriolis term is not taken into consideration as given 
in section 4.3.3. From the expressions for stability curves it 
can be seen that if the Coriolis term is neglected the stability 
characteristics of the blade tip motion will not depend on initial 

, . This is expected since the stability of linear 

systems should be independent of initial conditions. Also when 

= 0) relative to the spinning 
system the Coriolis acceleration does not effect the stability regions. 

An interesting effect of Coriolis acceleration is the creation 
of additional region of instability at a^ = 1/9. This region is 
not predicted by the linear analysis of section 3.3.2. 


there is no initial motion (Aj_ = 0, 



Note that A^^ 
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1/9 0.24 1/4 0.26 1 

DIMENSIONLESS FREQUENCY 

Fig. 5 Stability Chart for Initial Displacement '’^|=0.05 
and for Different Initial Speed 



1/9 


1 


0.24 1/4 0.'26 

DIMENSIONLESS FREQUENCY 

Fig. 6 Stability Chart for Initial Displacement 
and for Different Initial Speed 




Unstable regions at = 1/9 and a^ = 1 are quite narrow 
for practical initial values ('Qi ^ < O.oS") and widen for 

increased values of initial condition pareuneters r^, and 
On the other hand, the unstable region at a^ « 1/4 does not depend 
on the initial conditions. This region was predicted by linear 
analysis and is unaffected from the nonlinear Coiiolis term. 

Studying the three stability charts it might be concluded 
that fcr moderate initial values and low prec ssion rates 
the effect of Coriolis force on motion stability can be neglected for 
blades undergoing gyroscopic motion. 

So far the stability of tip motion has been analyzed by 
the perturbation method. Because of the uncertainty in the accuracy 
of this method it was decided to verify the stability charts presented 
in Fig. 4 , Fig. 5 , and Fig. 6 by the numerical procedure that was 
described in section 3.3.3. The solutions to the linearized equation 
(43) are stable if 

I V, (apTT) + Vj'(2pTr) I <a ( 63 ) 

where v^^ ( 6 ) and V 2 ( 6 ) are numerical solutions of equation (43) 
with the following initial conditions 

V, (o) =1 , v/(o) ■=■ O 

Vi(o)“ 0 , - O 

(64) 

Using this numerical method the transition curves that correspond 

to =0.1 are chec)ted. The stable and unstable 

points that are predicted by the numerical method are mar)ced on the 
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stability charts of Fig. 4, Fig. 5, and Fig. 6. They indicate 
that the stability curves obtained from the perturbation method 
are quite reliable. 

3.3.6 Forced Vibrations 

The equation of motion (15) indicates that the gyroscopic beam 
motion is a forced vibratory motion. From the stability analysis 
it has been determined that the motion in the unstable regions of 
the ^ ) plane is unbounded. This is due to the fact that 

the forcing function initiates vibrations which lead to parametric 
resonance. Hence it is only necessary to give the forced response 
away from the resonance regions. 

In the stable regions the solution to the equation of motion 
(15) can be assumed to be 

CK ^ Ul 6U (65) 


where Uj^ now includes the forced response 


*+ Pc Ul = e F “S.© -4- F Co'f p s ^ 


( 66 ) 


where F is given in equation (16) and is inthe order of 0(L), 
and u is the perturbed solution due to the remaining terms in (15) 

expecting the q^ term. As indicated before the terms with q^ has not been 
included in our analysis for the reasons noted earlier. 

In the stable regions the perturbed motion is bounded. 
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Fig. 7 Comparison of Numerical and Approximate Solutions in the 
Stable Region with Zero Initial Conditions 
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Fig. 8 Comparison of Numerical and Approximate Solutions in the 
Stable Region When Initial Conditions are u(0)=2.54 cm 

and ~= 2.54 cm 
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From equation (65) the error in taking u equal to u is in the order 
of 6 • Therefore in the stable region the forced response of the 

equation of motion (15) with error of order € can be obtained 
from equ^ttion (66) directly 


a - Uu * Ac«s\fTo^ •+ 

4- sio 2 0 -f coi’ & ^ 

Qo-'f Qo^l 


where 


A * u(o) 



vTs; 


rdu 1 26 F 

2 6 F Cot ^ 

^ L-c 

- \ 


(67) 


(68) 


The peak value of (67) can be determined by summing the 
amplitudes of individual harmonics: 


\ ^ “ \ /\ •+6 


do- H 


26 F cct ^ 

Qo - I 


(69) 


To check the validity of the solution (67) the equation of 
motion (15), disregarding the q^ term, is integrated directly. The 
predictor-corrector numerical scheme of Adams-Bash forth is used again. 
Both the numerical solution and the approximate solution are plotted 
in Figs. 7 and 8. Fig. 7 is constructed with zero initial conditions 


whereas in Fig. 17 the initial values are u(0) * 1 in.. 


q 

ae 


- I »n. 


e* c 
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gtu j 

which corresponds to c|4 It^o ' ^ in/sec. The solutions 

in Fig, 7 and Fig, 8 correspond to a stable point in the (a^, ^ ) 

plane which is a^ ■ 2 and 6 « 0,01, The particular blade data 

used in these figures is the same as those used in the stability 
chart. 

The approximate solution of (15) in the stable region which 
are obtained from equation (67) are in very good agreement with direct 
numerical solutions of (15) as can be seen in Figs, 7 and 8 , For 
smaller and more realistic 6 values the accuracy of the approximate 
solution will be even better. 

The previous discussion leads us to conclude that for practical 
purposes the pea)c value of the vibratory motion which is governed by 
equation (15) can be predicted from equation (69) provided the parti- 
cular blavle and rotary data correspond to points inside the stable 
region of the (a^, € ) plane. It is worthwhile to note that bending 

moments associated with pea)c tip displacements in this case ( € - O. ^l) 
are acceptable. Therefore the forced motion of a blade under gyroscopic 
disturbance is generally less important than the corresponding para- 
metric stability problem. 

3.4 Blades with Different Geometry 

In this part we show how to modify the analysis of previous 
section to take into account some common blade geometry, namely 
tapered cross sections and shrouding of blades. It can be shown 
that these properties do not change the general form of equation of 
motion (15) . The parameters in equation (16) are changed however. 
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They include blade geometry effects. Since the differential 
equation of motion (15) preserves its original form the stability 
analysis of previous section remains valid. 

3.4.1 Tapered Cantilever Blades 


In section 3.3.1 we have given eq :ation of motion for uniform 
blades that have a taper parameter, t^, « 1. If t^. ft 1 the parameters 
in equation (16) ta)(e the following form 


e = - Lot 
6 = JX/lo 

M * O.I%4 (ir - l) 

■ L 

« t (^r- 0 [TT^tr"^- 3 (fr “ l) ] 
7T^ 

-^0.0ih2 •+ ^.3bb^ ^ cos p 

Q. = 7^ 0?37 o.M3HceS^6 

N ( L 

i- (*fr ' 0 (- 0. o^q -f 0. IM3I o.oqa cos p )] 

j ~ — Sto 2 ^ (jnc^an<^ e<J ) 


] 


(70) 
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“ ^ 0.H3I -f ^ o.ndH Cos^^ 

4 (4r -0 (0.117 1 S- 0.1431 4 O OS3C0S^^ 


V- 


- -_L Sin ^ 

N L 


[o.aH«2 4(4r~0 


F--h£ll [o. 

K t ^ 


1343 


^ ^ 4 0.I%I7~ (70 cont.) 

(o.lofc4 4 0.1344“ ) ] 


The new blade parameters which appear in the above equation (70) 
are material density, f , the width of the blade, b which is assumed 
constant, the root cross sectional a^ea, A and t^, the area ratio of the 
tip to the root. N is some shape factor which is defined above. 

3.4.2 Tapered Blades with a Shroud at z = R + 3L/4 

Except the assumed displacement distribution (13) the general 
procedure described in Section 4.2 and Section 4.3 is exactly 
repeated. The assumed displacement shape is now 

5 ,4) ^ 4 (-f) ~ 3 j (71) 

which satisfies the geometric boundary conditions since 
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( 72 ) 


u(t) is again the tip displacement at z = R + L. 

For this case the parameters in equation (16) of cantilever 
blades take the following form 


6 » ~ 

e -0-/0) 

ISj = 0 0%57 (4r ~ 0 


Qo = 


ea 


^3.S33 (-fr-0^-^13*^ Ifc.S-(-fr-l) -I- vj 


^0,0571*) -t- O.O^S7~ ^ CaS ^ -{r (o.9^7f 






(73) 


^ ^00^71^ ~aN-+-Ncos^^ -^o.osrT-p 


u 


-t -fr ( O.^S-?! 


O.OSin 0 0 %^ 1 - + NCpS ^ 

-(r (O HS'?! } 

(continued) 
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(73 cont) 


3 . 5 Time-Varying Precession 

In Section 4.2 we derived the kinetic energy expression for 
space-fixed, constant magnitude precessional angular motion of the 
rotor. In practice, however, time dependent precessional velocities 
are very possible. In this section we will examine the case when 
the magnitude of the precessional angular velocity vector is an 
harmonic time function 

I J? I = cos (74) 


Under this condition the angular velocity vector is expressed as 

CO ' 4 _fLCoS>^ (coS^j 4 S(o9k) (75) 

where ^ is the frequency of the precession and 9= — oof' 
as before. 

3 .5.1 Line<-r Equation of Motion and the Its Stability 

The equation of motion is obtained by the same procedure that 
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was discussed in Section 4.2 and Section 4.3. Disregarding the 
nonlinear terms the equation of motion is given 


^ Qc + Z l 6^ -T ^ ~ [tos(l -F cos( i-’o)e3 
a© " a 2. «- 

9 

^ ^ - 0 -F ^ Co52(|-»-T)) ^ 

^ |ccs2(.-.»e]]u 
= - € F fan ? [ S(Y O +D) e Sio { \ - D') ©I 

: ir7e ( 1 +T)) e -F 

S|V< 2 ( ) - ■23') ^ ] 


where u indicates the tip displacement again. The quantities 

®o' ®i' ^1' ^2' ^ related to the blade parameters as follows 


0 ~ - tA)t 


e ^ 


J~L 

00 


<< I 


-0 ^ 




ET 


(77) 
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(77 cont) 


The equation of motion (76) represents a forced oscillation 
whose stability characteristics are determined by its homogeneous 
solutions. The homogeneous equation can be rewritten in the form 


• 4 ^ -hqoU=^u • ^ ( e , e) ^ 
dd 


(78) 


Assuming the solution of the above equation to be in the 

form 


U - A CoS (nTST a •+ o') ^ C = co<>st'Qnt 


(79) 


equation (78) may be solved by first substituting equation (79) 
into equation (78). The amplitude in equation (79) is than deter- 
mined by the equation 


1 iA 

A da 



e ( f.V)[(aN[To -F I -F T})e 


3 -F [(3xf^ — I 


~i3) a -F 


4 Stn[(' 2 vJ^ 4 I -r>)& -4- Sirs [C3>T^ - I -Fl))a ^c"]) 
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£ 

-*-2c] ->■ 5rrt[c2>T5; - 20)d -^;?c]) 


4 ^ 


fQo^ 


<^2 (s>n [(sJ^ 2) 0 7C~\ 4 sio f(2>Tq^ ~ a) © 


£ 

4 — ^ . - (fa l^(' 2 >)Te -*• 2 -*- ac] -4 5 io [( 2 Jao - 2 

- a.'?)') © -»- 2c^ 


sm [(2 >Iao -*“ a - 2-^)6 ac] *+ s.'o [(asT^ - a -»-a' 0)6 a cl) 


(80) 


If any term on the right hand side of this equation is not harmonic, 

the solution amplitude in (79) will be unbounded. This leads to 

unstable solutions at certain a values. The unstable values are 

o 

indicated schematically in Fig. 9. 


4-4 1 v-AA; 1 1 

0 T)^ 
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Fig. 9 Schematic Indication of Unstable Values 
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The general stability of equation (78) can be studied 
conveniently again by means of the parameter plane ^ ) • The 

plane can be divided into regions of stability and instability by the 
so-called boundary curves or transition curves, separating these 
regions. The perturbation method is used to derive these stability 
boundary curves at each of the unstable a^ values given above. 

In order to apply the perturbation technique for the determination 
of stability it is required that the function f( 0 , £ ) in equation 

(78) be periodic. 

The sufficient condition for f( 0 ,£ ) to be periodic is 

that 7 ) is a rational number, i.e., obtained by division of two 
mutually prime integers m^ and m 2 » 

(81) 

m, 

Under this condition the period of f ( ® ^ ) is 7T m^^ if m^^ and 

m 2 are both odd. Otherwise the period is 2 "TT m^^ . Later in this 
section we will restrict m^^ and m 2 to odd numbers only. 

By limiting f( 0,£ ) to periodic functions Floquet theory 

for linear differential equations with periodic coefficients is now 
applicable. Here we use the method of strained parameters again. 

The solution of (78) is assumed in the form 

u * Uo Ga, ■*■,.. (82) 

with the stipulation that the solution u is periodic with the period 
of f ( 0 , ^ ) , whereas the fundamental frequency a^ is given by 
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= bo ^ • . . 


(83) 


Substituting equations (82) and (83) into equation (78) a system 
of equations is obtained 


ao 


1 \doOo " 




de 


- b,Wo - ^ (Ac [co5(|-*-73)e -*■ CoS ( 1 - •o')© ^ 




de 


I bo(Jo “ Id, a, 4 b^Uo [coS (i -fo)© 4- CoS ( I - o) ©3 

d-'^CIo Co 5 27)0 <JaCoS2© 4 -^iJaCoSaCj -*•' 0 ) 9 

4-^ CoS ^ ( I - V)G 3 


(84) 


Equations (84) are solved recursively. Following the same 
pattern that was described in Section 3.3 secular terms are eliminated 
at each step and the equations for stability transition curves are 

obtained. We give the equations for these curves at Oo-lP’ , Qo“ I 
oi s ' j I -♦''0 

At q o = 

L “2 ^ J L?(i-u*) ^ J 

q, = r- iiii ^ \) “I r_±l__ ^ 1 

( 85 ) 
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Qc = \ 


La ^ 0-D^Xs--u') -* ^ J 

^ (i- (<1 •' ~ L ^(1“^’) ^ J 


3L 


€ 

(86) 


-D 


At q« = 




-t ^ 


(87) 


At 9, 


1 -+ *0 


H H 


( 88 ) 


Other transition curves are obtained accordingly. 
3.6 The Effect of Damping 


So far we have not considered damping in the analysis. In 
reality blade motion is damped by external and internal dissipative 
forces. Air resistance is a typical form of external damping and 
internal damping may be due to dissipative stresses in the blade 
which is called material damping. 

The dependence of dissipative forces on blade motion is 
quite complicated. In this section we shall study a simple model 
in which damping force is proportional to the velocity. Introducing 



ft damping force ^ to into the linear equation of motion 

(17) the equation becomes 

-^e-Q.C^S© -^€'<?aCo52e] a (89) 

die ^ d& '■ *• 

The constant c in the above equation represents the viscous damping 
factor. 

In order to study the stability of the motion governed by 
equation (89) the equation is reduced to a Hill's type of equation 
through the transformation 

U(b) e u(e) (90) 


under which 



(91) 


From the results of section 4.3.2 it can be inferred that the 

I 

solution to equation (91) is singular at Qo ^ * ‘A and 

Qo — ^ other wcrds, the singular points on the a^ 

axis of (a^, € ) plane for the u.^damped case are shifted to the 

left. In general c* is a very small quantity so that the singular 
points may be assumed approximately coincident with those of the 
undamped system. Since the other parameters in the differential 
equation (91) are the same as those of the undamped system of 
equation (17) the stability transition curves will have the same form. 
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In the unstable regions the solution to equation (91) has the 
following form 

Ji 6 

u($)^e $( 9 ) (92) 

where $ (6 07 t) - ^ , i.e. $ is a periodic function 

of ^ . The coefficient governs the growth rate of the solution 

u in the unstable regions and therefore sometimes it is termed the 
"negative damping coefficient". In the following steps we will 
obtain ^ in terms of the parameters in the differential equation 
(91). 

Inserting the assumed solution (92) into the differential 
equation (91) we obtain a differential equation in terms of the 
periodic function 

-V ^ -4. [Qo - € q, + 69 , Cose 4 € OxCos2e1 $ 

de^ ' die ^ ^ ^ 

= 0 

Following Whittaker's perturbation method that is given 
in reference ( 7 ) ws expand the solution and the parameters of 
equation (93) in terms of ^ as follows 

Lu 


(93) 


( 94 ) 
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Substituting the above expansion forms into equation (93) 
and rearranging in powers of € we obtain 

^ bo o 

be 1^1 = <^o - b, f2<e ~ <1, ^ 

0a -♦’bo^a •“ bi jZf/ 

- Co56^, ~ (JaCoSQ© 0O 


(95) 

Equations (95) should be solved in the neighborhood of the 
singular points mentioned before, name'.y Oq - and 

C(o — ' — I I • The resulting secular terms are eliminated 

(A» 

in each power of 6 so that periodicity of $ is maintained at 

^ I 

each step. At a. - * 'A we obtain the following 

conditions for eliminating secular terms in 



(96) 


(97) 


where A^, are the solution constants of ^o • Equation (96) 
represents a relation between b^^ and which can be used to 

express ^ in terms of the other parameters, A better approximation 
is obtained by eliminating secular terms in 0,, , This requires that 
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?-b: 




( 98 ) 


(99) 


from which it can be shown that 


O 


aoo) 


and 






(101) 


In equations (100), (101), and (96) the frequency expansion 
parameters are related to the negative damping parameters 
in the neighborhood of a^ * 1/4. We will give the expression for 
a constant ^ curve in the unstable region near a^ = 1/4. From 
equation (94) 

. ^ 

' ( 102 ) 

and 



Po - 



1L 

60 


J_ 

H 





(103) 
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Inserting equation (102) into eq. (101) and eq. (9fi) and then 
by substitution into equation (103) we get the expression for a 
constant curve in the (a^, € ) plane. 










(104) 


By finding the extrew*» of the above expression it can easily 
be shown that the lowest value of ^ is 


^ (105) 

' q, 

Thir equation can be used to determine the lowest 6 value 
that will cause instability when damping is present in the system. 

The analysis is repeated in the neighborhood of a^ = 1 to 
obtain the negative damping curves that correspond to a •* constant 

value. In the first power of 6 it is determined that 


/ 1 » 

and in the second power of 6 , 


(106) 


9' — • \/ 9,^ a. 

Equation (96) takes the following for.n 

T 

Po - = 1 -V € b, (108) 

LaJ 
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( 109 ) 


Inserting equation (109) into (107) and in turn substituting 
into equation (108) we obtain the negative damping curve that 
corresponds to ^ = constant curve in (a^, € ) plane in the 

neighborhood of a^ = 1 


qo ' 


cl ^ 




( 110 ) 


The minimum value of € on this curve is 


% 




(Ill) 


It should be noted that the negative damping coefficient in 
the neighborhood of a^ = 1 is in the order of 6 whereas that of 

a^ = 1/4 is in the order of € . Hence the growth rate is slower in 
the unstable region near a^ = 1. It can be expected that the growth 
rate at other singular points will be successively smaller. 

So far we have analyzed the growth rate of u that was defined 
by equation (90). The next step is to link the growth rate of u to 
that of u which is the solution to the equation of motion (89) for 
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blade vibrations whan damping is present. It must be noted that 
this dissipative damping has no relation to the negative damping that 
was mentioned before. 


Equations (92) and (90) can be combined to give 
(a{ 9) - e $C©) 


( 112 ) 


where $ is periodic. From this equation we clearly see the 
criteria for stability of damped blade vibrations. The vibrations 
with damping will grow in time if the n'egative damping coefficient 
associated with parametric excitation is larger than the dissipative 
damping coefficient of the system. 


c 

(aJ 


> 




(113) 


The stability transition curves of equation of motion (89) 
are obtained by setting 



(114) 


Inserting equation (114) into the negative damping curves 
(110) and (104) the stability transition curves for the differential 
equation of (89) are obtained. 


1 






CIS) 








il 

7 . 


(116) 
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These equations for transition curves could be obtained by 
directly applying Whittaker's method to equation (89) and search- 
ing for periodic solutions. 


Equations (115) and (116) are plotted in Fig. 10 for a 

particular blade data that yields = 1, q^^ = 1, and q 2 = 2. 

The stability transition curves that correspond to a particular 

V -3 

damping factor, oJ = 10 , are compared to the undamped case. 

If cO is 5000 RPM = 10 corresponds to a viscous damping 

factor, c = 0.5. As c becomes smaller the transition curves in Fig. 
10 approacn the curves for the undamped case. 


From Fig. 10 it can be seen that a finite damping coefficient 
causes blade instability to occur at lower processional rates. For a 
particular damping factor the range of safe precessional speeds can be 
determined from equations (111) and (105) when negative damping 
coefficient ^ is set equal to the damping factor, ^/oO 
Hence we obtain 






(117) 


and 


^ ? Cu3 

-J7.S ~ “ I 


(118) 


where Sis represents the safe precessional speed. From these 
equations and Fig. 10 it can be inferred that higher precessional 
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rates are required at a^ * 1 for instability as compared to 

a^ = 1/4. For the particular example (c = 0.5) the safe range 

of precession in terms of 6 - -^/ca> are = 0.002 

at a = 1/4 and Gs = 0.051 at a =1. 
o o 
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4. Arbitrary Blades Under Combined Rotor Spin and Space-Fixed Precession 


4.1 Introduction 

The stability analysis performed in section 3 are based upon 
an assumed blade deformation shape (see equation 3.13) which might 
be different from the actual shape that would occur. Another 
deficiency of the approach used in the last section is that pretwisted 
blades can not be analyzed in the same fashion. In this section a 
general numerical formulation to analyze the stability of a pretwisted 
blade of arbitrary cross-section is presented. Results obtained 
by this method are compared with those computed in section 3. In 
addition the behavior of a turbomachinery blade with realistic 
dimensions is studied. The same approach may be used to analyze 
forced vibrations of a pretwisted blade under gyroscopic excitation. 

The first step is to derive a stiffness matrix for the 
twisted blade. In order to handle blades wivn arbitrary geometry, 
finite element approach is selected. For a twisted blade usually 
the shell element is the most appropriate finite element. However, 
the numerical method for dealing with the stability problem is 
quite complicated. It is desirable not to employ too many degrees 
of freedom for the modeling of the blade. For this reason the beam 
element is chosen. Since the problem interested in this work happens 
almost exclusively to long and slender blades this simplification 
is not considered serious. A 12 x 12 beam element stiffness matrix 
for twisted blades is formulated based on a theory proposed by Downs 
(10) . Static deflection solutions are used to study the accuracy 
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of this stiffness matrix. 


The derivations for mass matrix and the changes in stiffness 
matrix due to rotational motions are based on Lagrange's formulation 
of equations of motion. Before the equations of motion with time 
dependent coefficients are analyzed the natural frequencies in both 
rotating and i on-rotating force fields are calculated. The computed 
solutions are compared with known analytical and experimental results. 
The subsequent stability analysis of the pretwisted blade under 
both rotor spin and space-fixed precession is studied via numerical 
methods based on Floquet's theory. 

4.2 Stiffness Matrix 

The elementary beam theory cannot be directly applied to 
the study of twisted blade. The modified beam theory assumes that 
the undeformed blade is pretwisted along its length about a straight 
longitudinal axis through shear centers of cross sections. These cross 
sections are parallel to the rotor fixed xy plane (see Fig. 3.1). 

The blade surface is supposed to be composed of helical fibers. If 
the position of an arbitrary point on he cross section is indicated 
by its principal coordinates ^ and , the distance of this point 

from the shear center is 

a a 

d- (r. ^ n) a) 


By designating oi as the initial twist of the blade in radians 
per unit length and considering a portion of the blade with unit 
length, any helical fiber of length j?n which is at a distance d 
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from the shear center forms a small spiral angle ^ (see Fig. 1) with 
the longitudinal axis. The length of this helical fiber can be 
approximated as 

|„ « ( I (2) 


Knowing this length, the axial strain of a helical fiber due 
to bending, axial load and torsion can be derived. The axial fila- 
ment strain caused by an extension of "a" per unit length is 



(3) 


Neglecting higher order terms in ^ this equation may be reduced to 



(4) 


The axial strain due to an elastic twist p per unit length can 
be given in a similar manner 


= 




Again ignoring higher order terms in ^ and 


(5) 




( 6 ) 
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The total axial strain form axial extension, bending and torsion 
can be derived 




( 7 ) 


where ^ ®x**^ ~ *1" ®y positions of the centroid)and 

where and Ujf are the deformations in the principal directions. 

The stiffness matrix is obtained from the potential energy 
of a beam element 




( 8 ) 


where ^ is the deformation strain and its corresponding stress. 

Substituting the axial strain given by equation (7) into the strain 
energy expression (8) the stiffness matrix is determined after 
performing the integration. Since torsional strain was not included 
in the strain equation (7), torsional rigidity has to be added to 
the final stiffness matrix. It is worth noting that other effects such 
as shear deformation, warping and rotary initial terms can be included 
in the strain expression (7) if they are considered important. The 
integrand in equation (8) may be written in matrix form 






--fC «<£>=■ ~ 
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( 9 ) 




For the finite element formulation the deformations within 

an element have to be expressed in terms of nodal displacements. 

This is accomplished by assuming reasonable distributions for the 
element deformations between end nodes and relating these deformations 
to nodal displacements by shape functions. Element bending deforma- 
tions and in equation (9) are in the principal directions 

of the cross sections and hence they rotate with the pretwist angle 
along the longitudinal axis. This will lead to complicated finite 
element shape functions. In order to keep the shape functions 
simple the principal displacements ZT and V at the middle 
of the blade element are selected to represent the element bending 
deformations. They are related to bending displacements of an 
arbitrary cross section by (see Fig, 3) 


6(£* U coSo^S ■»“ N/sineiS 
Urf ^ -ZTS'Oc^S V/CoS^S 


(10) 


The coordinate s originates from the middle of an element. It 
coincides with the longitudinal z axis and they have the following 
relationship 

0 0 

Sr B - P - 2 ■^'5 (1 

where n is the number of the element and reprefents the length 


of each element. The length of the nth element under consideration 
is J?. R is the rotor radius. 



All the displacements are approximated by polynomials 


(X) * C,S Ca 

M * S Cs S “X C(, 5 

The constants are determined by the nodal displacement variables 
at the two end nodes (1 and 2) 


f ^ u). j 0*^ 


(13) 


The expressions of are listed in Appendix A. 

The deformation variables used in equation (9) can now 

be expressed in terms of nodal displacement variables . '^hfi 

bending curvatures in ^.P^VL^ are related to curvatures in the 
middle of the element. 


(14) 


ci O-'^ • j X d \j ic 


Differentiating equation (12) with respect to s the following 
equations are obtained 
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Cm +aCsS + 3C*S* 
* ac, + tCiS 
^ = C« + aC,s 4- -aCioS* 


( 15 ) 


The constants are known in terins of nodal displacement variables 
M • They are substituted into equation (15) , and with equation 
(14) the transformation matrix between deformation variables 
and nodal displacement variables is obtained 

\ lAI. 

dl 

dz 





(16) 




or 



J 


are given in 


The elements of transformation matrix 
Appendix B. 



Inserting equation (16) into equation (9) and after performing 
the matrix algebra the potential energy is 



(17) 


where the matrix [S] is 


[sv- 


(18) 


From the potential energy equation (17) the stiffness matrix is 
derived 

“'v 

To extend this integration for the stiffness matrix to a 
relatively complicated cross section of the type which is liable 
•■;r occur in real turbomachinery blades requires the introduction 
of a method )cnown as the isoparametric transformation in the finite 
element analysis (^) . With this approach the evaluation of the 
stiffness matrix represented by equation (19) can include cross 
sections of any arbitrary form. This involves the distortion of 
a simple basic square form into areas with arbitrary boundary shapes. 
The accuracy of the 'mapping* depends on the number of nodes 
introduced around the circumference. In this report all results 
are based on eight boundary nodes. 



In the standard isoparametric formulation of three dimensional 
distortion, mapping takes place in all three directions. In the 
present approach, the integral (19) of the stiffness matrix already 
contain the blade pretwist effect. Therefore, it suffices to perform 
distortion in two dimensions of the blade cross-sectional area only. 
This reduces the amount of required numerical computation in evaluat- 
ing the stiffness matrix. It also assumes that the axis of pretwist 
is a straight line. 

The isoparametric formulation begins with the establishment 
of a one-to-one correspondence between the Cartesian and curvilinear 
coordinates. The variables £ and '(( in matrix [P] of equation 
(18) are related to the curvilinear coordinates x and y 



Once such coordinate relationships are known, shape functions can 
be specified in local undistorted coordinates and by suitable 
transformations the stiffness matrix evaluated. 

The most convenient way of establishing the coordinates 
relationship in (20) is to use the shape functions associated 
with the eight boundary nodes. This may be written as 

g - Mi§, + NaSj . Ntit 

M,<}, .... 


( 21 ) 



in which are shape functions given in terms of the local 
coordinates (x,y) . Based on the basic square form the shape function 
expressions are 
Corner nodes 


Ni’-^ +tjljj - l) (22) 

Mid-side nodes 

X ;-0 , N,- + ) 

iji^o , Ni ~ 'j ' 

where i is the number of the shape function, and x^, are either 
+1 or -1 or zero depending on the location of the boundary node. 

To evaluate the stiffness matrix in equation (19) two 
transformations are necessaiy. In the first place C and are 
defined in terms of local (curvilinear) coordinates x, y by equation 
(21) . In the second place the volume over which the '.ntegration 
has to be carried out needs to be expressed in terms of the local 
coordinates. A standard process will be used which involves the 
determinant of the Jacobian matrix [J] . Thus 


^ cJet 


( 24 ) 


The Jacobian matrix [J] is found explicitly in terms of the 
local coordinates by 



( 25 ) 





' 5, T7, 

) 

I 

€, n% 


After substituting equations (21) and (24) into the integral 
express ‘on for the stiffness matrix (19) the stiffness matrix of an 
arbitrary cross-section is 


[s ] ci«t 


(26) 


The above integration must be evaluated numerically. A flow-chart 
showing how to compute this 12 x 12 matrix using Gaussian integration 
is given in Appendix C. 

In order to establish the accuracy of this numerical method 
the static deflections of cantilever beams with arbitrary cross- 
section were investigated. For different elliptical and triangular 
cross-sections the errors in computed free and deflections using 
this approach with six elements remained below 3% of the theoretical 
values. This seems adequate for most engineering applications. 

In particular for the stability analysis of interest in this wor)c 
the accuracy of the approximated stiffness matrix is sufficient. 

Another interesting example selected to evaluate the numerical 
accuracy in computing the stiffness matrix is the cross-section 
defined by a convex and a concave edge and shaped like a crescent 
(see Fig. 2). Two nodal distribution schemes were tested. They 



are shown in Fig. 2a and 2b. For a 6 inch long cantilever beam 
with 3000 lb. load at the free end the vertical and horizontal 
free end deflections are compared with those obtained by the 
numerical method (six elements) : 



1 reckon 


A 





. 0 %^% 

Bea/vi SoU'boo 


. 0%11 


This indicates that the eight nodes in case b are better placed 
than the ones in case a. The numerical interpolation is quite 
sensitive. If three nodal points are employed along the circular 
edge only, the results can be as much as 100% off the correct values. 
This is due to the fact that the three node parabolic approximation 
cannot represent the circular edge correctly. 

For twisted blade with rectangular cross section some 
experimental results were presented by W. Carnegie (^) . In 
these experiments the blades were uniformly pretwisted over a 6 inch 
(15.2 cm) length and had cross sectional dimensions of 1 in x 1/16 in 
(2.54 cm X 0.16 cm). The blade material was mild steel with density 
of ^ = .000735 lbf“sec^/in^ and with Young's modulus of E = 30 x 10 
psi. Using the same blade data and six finite elements the static 
deflection under unit load at the tip was computed. The results 
are compared to Carnegie's experimental results in Fig. 4. As can 
be seen deflections in the horizontal and vertical directions become 










due to pretwist. Satisfactory agreement is obtained in these 
comparisons. 


4,3 Kinetic Energy and Mass Matrix 


The mass matrix is der ‘ j ^ " rom the kinetic energy which 
also determines the changes in the stiffness matrix due to rotational 
motions. For the case .. % space fixed precession these changes in 

stiffness matrix are time dependent and might lead to instability 
of the system. Without precessional angular velocity usually the 
primary interest is in the calculation of natural frequencies. 

After deriving the kinetic energy and the associated matrices, the 
correctness of these matrices will be tested by comparing the 
computed natural frequencies in rotating and non-rotating force 
fields v‘th known experimental and analytical solutions. 


Beginning with the displacement vector P of an arbitrary 
point (x,y,z) of the beam 



( 27 ) 


The x,y,z coordinates employed in this displacement equation are 
the rotor- fixed coordinates shown in Fig. 3.1. The corresponding 
velocity is 


■ 5 = 

- U\^^) J If ■ ^ Js) k + W X r 


; r *'iv >'r> 




where 6J is the angular velocity of the rotor and was defined in 
equation (3.1) of section 3. The velocity V is to be inserted 
into the Jcinetic energy 


T = i? I v-v 


For the finite element analysis, the deformation variables 

c C ^ OV *5 


in the velocity V equation (28) have to be expressed in terms 
of nodal displacement variables employed in the stiffness 

matrix (19) . Due to the angle between the major principal bending 


axis of a cross section and the rotor spin axis the deformation 
variables can be related to deformation variables at 



( 30 ) 
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or 


\ - [I2] 


where the angle ^ is between the rotor spin axis and the major 
principal bending axis of an arbitrary cross section. If the pretwist 
angle per unit length ^ is constant the angle ^ is (see Fig. fi) 




(31) 


with as the blade setting angle at the root. 

To be consistent with the shape fxinctions employed in the 

— A 

derivation of the stiffness matrix the principal displacements U 
and \/ at the middle of an element are selected to represent the 
element bending deformations. The angle of twist of the nth element 
middle cross section is 
n-» 

pc (32) 

J *1 

where o^i and are the pretwist rate and length of elements 

before the nth element. With in equation (30) the 

deformation variables ^ 6 (f^ functions of 


,W.y 


' 5 ? 


^ 0?? 


( 33 ) 


Repeating the procedure used to derive equation (16) for the 
stiffness matrix the above variables can be determined in terms 
of nodal variables 
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( 34 ) 


The elements of matrix fr*l are listed in Appendix D. With 
) in equation (30) and from equation (34) 

the variables in the velocity V equation (28) are expressed as 
functions of the nodal variables . Substituting this velocity 

into the )cinetic energy (29) , this leads to the Jcinetic energy of 
an element 

*+ ^ M4^ ^ Jv (35) 


where are the nodal displacement variables and their 

time derivatives. The matrices [Mil are 

[M,l= [T,l7T,y[A.l[T,][T,] 

[M.>[T,r[T,f[A.][T,llT,] 

[Mil *hT[T,r[A,] 

?M.V[T,nT.fjA.KT.lW 


The matrices 


are listed in Appendix E. 



with both kinetic and potential energies the equations of motion 
can be derived via Lagrange’s formulation. The mass matrix is obtained 
after integrating matrix fM.T . This mass matrix is not exactly 
the same as the standard mass matrix because of the pretwist effect. 

The differences are, howe^ er, very minor. The elements of matrix 
[mO are time dependent. They reflect the effects of both 
rotational and precessional motions on the beam stiffness. This matrix 
is combined with the stiffness matrix of equation (26) to form a 
time dependent matrix representing the spring stiffness coefficient 
matrix of the equations of motion. Matrix does not appear 

in the equations of motion. The forcing terms of the equations of 
motion are obtained from matrices and The integrations 

in equation (35) are performed numerically. For arbitrary cross-sections 
the same numerical distortion procedure used in the evaluation of 
stiffness matrix has to be applied. 

Without the precessional angular velocity the integrated 
matrix tM.l represents the change of blade stiffness due to 
rotor spin. To demonstrate the accuracy of the present approach 
the following numerical exeimple was calculated; length, = 1 in., 
width, w = 3 in., thickness, t = 0.09 in.. Young's modulus, 

E = 1 x 10 psi, Poisson's ratio, *0 = 0.3, and mass density, 

^ = 2.587 x lO”^ Ibm/in^. The results are presented in Figs. 

5 and 6 and are compared with the experimental and numerical results 
obtained by MacBain (J^) using NASTRAN program and 230 plate elements. 

The present numerical results, based upon the beeun theory, are 
obtained using 6 beam elements. In Fig. 6, the first, second and 
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Nondimensional Frequencies, ^ 
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Rotational speed 

Fig. 7 Natural frequency vs rotational speed 
(30 deg pretwist angle) 



third bending frequencies and the first and second torsional 
frequencies at zero spin velocity are plotted vs. the total pretwist 
angle of the blade. In Fig. 7 , the same frequencies with 30° deg 
pre twist angle are plotted vs. the rotational speed. In this case 
the blade setting angle = 0 and .ub radius R = 7 in. 

The same examples were used by Chen and Dugundji (14) to show 
the accuracy of a twisted beam finite element derived from the 
governing differential equations of motion provided by Houbolt and 
Brooks (^) , Since the theoretical assumptions made by Houbolt 
and Brooks are essentially the same as those made in this report 
there is practica"* y no difference in the numerical solutions. 

Comparing the derivation of stiffness and mass matrices for the 
finite element analysis the present approach is superior because 
of i'' s straight-forward simplicity. Nonlinear material r geometrical 
effects can be easily included in the derivation. The evaluation 
of cross sectional properties are based on the shape of the area and not 
its moment of inertia which is often difficult to determine. 

4.4 Stability Analysis 

The homogenous equa - iS of motion for the study of stability 
due to precessional motion can be formulated as a system of linear 
differential equi^tions with periodic coefficients. The system 
has the following form 

^ [K*l \u,\ = o (37) 
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is the mass matrix defined by equation (36). [K*J is the 
effective stiffness matrix which is the combination of the structural 
stiffness matrix [K] defined in equation (26) and the [M 2 ] matrix 
given by equation (36) . The periodicity enters equation (37) 
through the matrix, [M 2 ] » whose time dependent terms are due to 
the precessional motion. 

A vast array of literature exists on the subject of ordinary 
differential equations w:.th periodic coefficients. Some of them 
are procedures to study the stability problem. After reviewing these 
procedures, it was determined that most of them are not suitable 
for systems with many degrees of freedom such as equation (37). 

For large systems a general method of solution proposed by 
Friedman ( 16 ) was found to be most appropriate. The details of 
the method are not presented in this report because the procedure 
in reference (]^6) was adapted here without any major revision. 
Following the suggested procedure and employing the improved 
numerical integration scheme based on the Runge-Kutta method a 
computer program was developed to study the stability regions of 
the equations of motion governed by the matrix equation (37). 

The first problem studied was the straight beam with constant 
cross section. This case was already analyzed using another method 
in section 3.3.2. The results were plotted in a stability chart 
in Fig. 3 of section 3. The linearized finite element solucions 
based on equation (37) with one and two elements were able to 
suplicpte the same stability chart. This demonstrates that the 
basic approach employed in the finite element solution is correct. 


Furthermore, since the stability chart is now obtained by two 
different solution methods, there is much more confidence in the 
validity of this chart. 

The next case studied was the pretwisted beam with constant 
cross section. The stability chart for this case is basically very 
similar to the stability chart of the beam without pretwist. The 
effect of beam pretwist is to increase the structural stiffness 
in the effective stiffness matrix [K*] of equation (37). This 
leads to an increase of the parameter a^ of the stability chart 
in Fig. 3 of section 3 depending on the amount of pretwist. Since 
the regions of instability decrease with increasing values of a^, 
the pretwisted beam is less prone to unstable dynamic motion under 
precessional rotations. 

The computation for each point in the stability chart using 
the finite element solution is very expensive and time consuming. 
Not enough stable cr unstable points were calculated to chart 
a complete stability diagram. 


5. Forced Vibration and Flutter 


5.1 Introduction 

The complete time dependent problem under precessional 
rotation solved by the finite element method can be reduced to a 
system of ordinary differential equations of the characteristic form 

(38) 


in which [M^ ] and (K*J are assembled mass and effective stiffness 
matrices. The /orce matrix {F}is obtained from matrices {M^} and 
defined in equation (36) . In general, the above equations 
are non-linear; only linear cases will be studied. The forcing 
matrix {f} contains functions of first and second order in spin 
velocity; hence it is periodic. In the following the forced vibra- 
tion problem will be solved by the determination of periodic responses 
to equation (38) . 

The availability of the twisted beam finite element 
program and a subroutine for unsteady subsonic aerodynamics (17) 
made possible a short study of "fan flutter" under realistic 
conditions of structural coupling due to twist. The effect of 
centrifugal forces on the flutter mode is taJcen into account as 
well. The terms associated with precessional rotation were removed. 
The equations of motion to study flutter have the following form 

[kM = o 

( 39 ) 



Matrix [C] represents aerodynamic damping and contains complex 
terms. Substituting 




( 40 ) 


into equation (39) the characteristic equation is obtained 


o 

where oC and {u} are complex. The real part of the solution 
represents a decaying vibration. The complex eigenvalue problem 
involved in Eq. (41) is solved by a numerical procedure given in 
reference (^) . 


5.2 Forced Vibration 


In section 3.3.6 the forced vibration of the simplified 
single degree of freedom system was studied. It was shown that 
dynamic responses are not sensitive to time dependent quantities 
in the stiffness term. It seems reasonable to remove time dependent 
quantities from the effective stiifness matrix [K*] and reduce 
the coefficients of the differential equations (38) to constants. 

The forcing term in Eq. (38) is periodic. It can be 
expressed as 

where o( is complex. A general solution can be given as 


( 43 ) 



Substituting (43) in Eq. (38) gives 

U * w) 1,5^ . Tf’i 

By inverting the matrix 

^ * MS) 

and pre-multiplying with the force amplitude {f} the dynamic response 
{u} can be determined. 

The first computational example selected for comparison is 
the straight beam without pretwist. Using a model with one finite 
element the same tip displacement responses were reproduced as those 
given in Fig. 7 and 8 of section 3. A few pretwisted beams were 
studied as well. In general, the peak displacements are smaller for 
pretwisted beam. 

5.3 Flutter 

The computational model chosen for examination is a steel 
cantilever rotor blade uniformly twisted from root to tip such 
that the stagger angle varies from 30 degrees at the attachment 
to 60 degrees at the free end. The blade is untapered and the 
cross section is taken to be a thin rectangle of 5.08 cm chord 
and 4% thickness ratio. Pitch/chord ratio at midradius is 0.878. 

The hub/tip radius ratio is 0.5. Axial inflow is assumed, uniform 
radially, such that the relative flow angle at the blade tip is 
60 degrees. This specification determines the air velocity triangle 
at every radius; the aerodynamics are computed at the mean radius 



of each beam (blade) element using the formulation outlined in [17] . 


The method of solution of the equations of motion was 
essentially the same as that described in reference [18 ] , and known 
as the ”p-k method”. For a given Mach number, a value of reduced 
frequency, k, is assumed, the unsteady aerodynamic forces are computed 
using the aforementioned subroutine. The equations are then solved 
for the complex eigenvalues, co = w + i cu . The motion having 
been assumed to have a time dependence described by exp(io)t), the 
computed value of w., provided an improved estimate of the reduced 
frequency, k = Uj^C/(2V) and the procedure can be iterated until o)j^ 
does not change. The converged value of gives an estimate of 
the flutter frequency; the associated value of at convergence 
is an indication of the nonaerodynamic forms of damping (positive 
or negative) that must be supplied to maintain the constant amplitude 
"flutter" condition. In the absence of such external damping, 
the flutter condition occurs with co^. = 0; thus Wj. < 0 indicates 
instability, or potentially divergent motion. 

In the previous series of calculations the interblade phase 
angle is a parameter. In practice this is limited to a finite 
number of possible values 2 tt n/N where N is the number of blades 
in the row and n is any integer, 0< n 5 . N. For the extraction 
of the greatest amount of information it is convenient to let the 
interblade phase angle range continuously over 0, 27r. Thus, at 
each Mach number, a closed contour is obtained in o)j coordinates 
with interblade phase angle as a parameter. 

The true flutter condition occurs at that Mach number where 
a closed loop is obtained, lying entirely in the positive Wj region 


expect for one point at which the contour ' tangent to the axis. 

Figures 4-1, 4-2, and 4-3 show the contours obtained for 
the model blade at absolute approach Mach numbers of 0.75, 0.85 and 
0.95. It appears from these figures that flutter will occur at a 
Mach number just below 0.85 with an interblade phase angle in 
the neighborhood of 100°. 

In each of the figures a cusp-like tangency or near tangency 
of the contour occurs in the neighborhood of = 1023 to 1253 rad/sec, 
depending on the relative Mach number, and hence rotor speed. 

These points, and the irregular nature of the curves near these 
tangencies, are attributable to the "aerodynamic resonance" 
phenomenon. At the combination of parameters which produce 
aerodynamic resonance the unsteady aerodynamic loading on the 
blades vanish [18] and the blades vibrate as if in a vacuum. Thus, 
the computed eigenvalue is real and corresponds to the rotating 
natural frequency ^ vacuo . Although the aerodynamic damping does 
not go negative, operation at this point could result in large 
response to forcing and may be shown in some instances to lead 
to the accumulation of fatigue damage. 

True flutter of the coalescent type is illustrated in Fig. 4-3. 
Here a slight amount of mechanical damping would produce a flutter 
point at (iJ_ = 2330 and 100°. The mode shape represents a 
combination of the modes associated with the two lowest eigenfrequen- 
cies at o)„ =1138 and oo = 3296. 

H x\ 

One concludes from this brief study, and other supplementary 
data, that the structural coupling introduced by blade twist results, 
at least in these fairly repiesenative cases, in a classical 
"coalescence- type" flutter. 
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ig. 4-1 Lowest Mode Eigenvalues for 0.75 Approach Mach Number 






Fig. 4-2 Eigenvalues for Coalescing of Two 




coo 


ORIGINAL PAG£ tS 
OF POOR QUALITY 


r 

o 

o 

' • 

fM 



Fig. 4-3 Lowest. Eigenvalues for 0.95 Approach Mach Number 


6. Conclusions 


The 3 year program of research supported by NASA and reported 
herein has resulted in four papers listed as references 1,2,3 and 
5. Additionally, an interim technical report listed as reference 
16, was issued after completion of the first half of the contract 
period and has been incorporated into the present report. 

Analytical representations coupled with computer-based 
numerical solutions have been completed for a number of significant 
blade dynamics problems. The earliest of these are point mass 
and uniform cantilever models. A beam-type finite element, dis- 
playing the important foreshortening property, has been developed 
and forms the structural basis for all of the later prograras. 

Blade twist, taper and cross-sectional shape are accounted for. 

• Stability of a blade on a spinning and precessing rotor has 
been studied using this beam element. Both steady (ramp function) 
precessional displacement and harmonic precession have been 
studied and instability has been shown to be possible in both 
operational modes. The critical parameters discriminating unstable 
motion are found to be dependent primarily on the ratios of 
rotating natural frequency to spin frequency and of constant* 
precession rate to spin frequency. The blade setting angle and ratio 
of blade length to attachment radius are also of significance. 

With harmonic precession the critical values of spin rate are 
found tc be dependent upon the frequency of the harmonic variation 
of the precessional motion. 


* In the case of harmonic precession the maximum rate (or amplitude 
thp constant of the steadv orecession case. 



• Forced vibrations, outside the regions of instability and 
excited at three basic freguencier (a subharmonic, the spin 
frequency and twice the spin frequency) are attributable to 
the interaction of precession and spin and may be calculated 
in a straightforward manner. 

* Blades (beams) modeled as pretwisti^d beam-type finite elements 
were subjected to su?joonic cascade flutter analysis. It was 
demonstrated that it is important to ta)ce into account the 
structural coupling due to pretwist and to properly model the 
effect of the rotational field on the rotor blade mode shape. 

The general conclusions of greatest importance are: 

1) the rotor blade model for dynamic analysis in complex 
rotational fields and in air aerodynamic environment 
should properly reflect foi'eshortening and pretwist 
for reliable predictive capability, and 

2) thin blades, on the order of 5% tip th’clcness, should 
properly be analyzed using a plate-type finite element 
incorporating foreshortening (as yet unformulated) and 
with arbitrary orientation in the force field due to 
arbitrary rotationaJ motion (spin plus precession). 

3) Twisted cantilever rotor blades at subsonic speeds tend 

to flutter in the coalescent mode in which the two gravest 
mode couple to produce the negative aerodynamic damping. 
This result is interesting since it is at odds with some of 
the theoretical predictions for untwisted blades. 
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